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Abstract 

The  most  common  solution  technique  is  subspace  iteration,  a  combina¬ 
tion  of  inverse  iteration  and  the  Rayleign-Ritz  procedure.  Some  diffi¬ 
culties  of  this  method  are  mentioned  as  well  as  ways  to  avoid  them  by 
means  of  spectral  transformations.  This  permits  use  of  the  Lanczos  algo¬ 
rithm  and  yields  significant  reductions  in  cost. 
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1.  INTRODUCTION 

In  order  to  avoid  excessive  generality  we  assume  that  K  and  M  are 
real,  symmetric  n  x  n  matrices.  The  goal  is  to  compute  all  the  eigen¬ 
values  lying  in  a  given  interval  [t  t^]  together  with  the  correspond¬ 
ing  eigenvectors. 

For  n  >  2  the  eigenvalues  will  he  real  provided  that  <K  +  pM  is 
positive  definite  for  some  choice  of  <  and  p  .  For  simplicity  assume 
that  either  K  or  M  is  positive  definite. 


Example.  If  K  =  M  =  | 
eigenvalue  belonging  to 


then  every  number  (real  or  complex)  is  an 


The  pair  (K,  M)  is  often  called  a  pencil.  Its  eigenvalues  and 
eigenvectors  are  denoted  by 


A  <  A  <  <  A 

1  2  •  ’ '  n  ’ 


Z1 »  z2 ’  * • • »  zn  • 


Techniques  for  solving  the  problem  when  n  <  100  are  described  in  [Parlett, 
1980]  and  will  not  be  mentioned  here. 

The  most  important  factor  in  selecting  a  particular  algorithm  is  the 
relative  cost  of  solving  (K  -  oM)u  =  v  for  u  to  forming  the  product 
u  =  (K  -  aM)v  .  Observe  that  this  is  simply  the  divide:  multiply  ratio 
for  the  given  matrix  pencil.  Linear  equations  can  be  solved  iteratively 
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(using  conjugate  gradients  or  Lanczos)  when  triangular  factorization  is 
not  possible.  The  remaining  sections  show  how  this  divide:  multiply  ratio 
affects  the  solution  technique. 

With  apologies  to  structural  engineers  we  use  capital  letters  for 
matrices,  small  roman  letters  for  column  vectors,  and  small  greek  letters 
for  numbers. 


2 .  SUBSPACE  ITERATION 

See  [Bathe  and  Wilson,  19?6]  or  [Parlett,  1980]  for  more  details. 

It  is  assumed  that  M  is  positive  definite. 

Before  the  start,  a  shift  a  is  chosen  in  [x^,  t^]  and  the  triangular 
factorization  K  -  aM  =  LDL*  is  computed.  Note  that  for  each  i  =  1,  ...,  n  , 


(K  -  CM)  "Slz.  =  z.  v.  , 


l  l 


Vi  ~  IxTiaT  ' 


[2-1 ) 


So  the  power  method  with  the  operator  (K  -  aM)  "Si  would  converge  to  the 
eigenvector  z^  belonging  to  the  closest  to  a  .  The  usual  normalization 


z.*Mz  =  5  (Kronecker 
^  <J  i  j 


Of  course,  K  -  aM  is  not  to  be  inverted 

explicitly. 

The  number  of  eigenvalues  in  [x^,  is  unknown  but  subspace  iter¬ 

ation  must  begin  with  the  difficult  choice  of  the  dimension  p  and  the 
somewhat  easier  choice  of  starting  vectors  (x^,  x^,  . 


x  )  =  X0  with 
P  0 


Xq*MXq  =  1^  .  The  usual  implementation  is:  for  k  =  1,  2,  ...  until 
convergence  repeat 
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1. 

2. 


Solve  (K  -  aM)Yk  =  MX^  for  Yk  . 

Compute  the  reduced  matrices  (or  projections) 


\  -  V°k  • 


\ '  Yk*mk  • 

3.  Solve  the  small  problem 

normalized  eigenvectors  g^,  g^  . 

'*■  s«  \  *  V*  • 

5.  Test  for  convergence. 

To  test  for  convergence  one  should  use 


=  0  for  eigenvalues  9^  and 
Set  C>k  ( gj. ,  •  •  •  ,  gp )  • 


Theorem.  For  any  x  ¥= 0  and  any  scalar  0  there  exists  an  eigenvalue 
of  (K,  M)  such  that 

|A  -  0j  <  ll(K  -  6M)xl  ,  / IlMxll  . 

M 


Here  lluiL  =  /u*Hu  .  Thus  IlMxll  ,  =  /x*Mx  which  is  easy  to  compute  but  the 
H  M 

cost  of  the  numerator  discourages  the  use  of  this  theorem.  Often  users 
wait  until  the  values  of  9  "settle  down"  to  the  desired  accuracy  but 
such  a  criterion  is  not  completely  reliable. 

A  serious  difficulty  with  subspace  iteration  occurs  when  (K,  M)  has 
eigenvalues  outside  both  ends  of  [ ,  t^]  .  To  illustrate  what  happens 
consider  an  extreme  case:  suppose  that  at  some  step  k  one  of  the  Ritz 
vectors  (the  columns  of  Xk  )  is  x  and  satisfies  x  =  (z^  +  z^)/^2  and 


A  <  T  <  t  < 
112 


30 


A 

1 


T2  * 


It  turns  out  that  the  corresponding  Ritz  value  is 
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Z1*M21,V23Q*MZ3Q,X3Q 

z^^Mz^+z^^Mz^ 


VX30  _  VT2 


Although  9  is  in  the  middle  of  [t  ,  T2^  does  not  signal  that  there 
is  an  eigenvalue  A  nearby.  R.L.  Taylor  calls  these  "ghost"  values  and  in 
[Scott,  198l]  it  is  shown  how  their  presence  can  prevent  the  discovery  of 
good  approximate  eigenvectors.  One  remedy  is  to  compute  all  eigenvalues  and 
eigenvectors  less  than  and  then  keep  orthogonal  to  these  eigen¬ 

vectors.  That  is  expensive. 


3.  TRANSFORMING  THE  SPECTRUM 

Steps  2,  3,  and  U  of  subspace  iteration  carry  out  the  Raleigh-Ritz 
procedure  which  delivers  the  best  approximations  to  eigenvalues  of  (K,  M) 
from  the  space  spanned  by  .  To  remove  the  troubles  of  Section  2  apply 
the  Rayleigh-Ritz  procedure  to  ((K  -  aM)  I)  because  the  wanted  eigen¬ 
values  v.  ,  see  (2-l),  are  the  extreme  ones.  However  (K  -  oM)  is 
1 

not  symmetric.  There  are  two  ways  out  of  this  difficulty. 

I.  The  Spectral  Transformation  (see  [Ericsson  and  Ruhe,  1980]  for  details). 
Factor  M  =  LL*  and  use  L*(K  -  aM)~^L  instead  of  (K  -  aM)  M  . 

In  other  words  change  (K  -  AM)z  =  0  into 

(L*z)  =  L*(K  -  oM)-1L( L*z )  .  (3-1) 


The  powerful  Lanczos  algorithm  can  be  used  with  the  new  operator.  Another 
advantage  is  that  the  error  bound  of  Section  2  tec  cries 
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I  6  -  v  I  <  “L*(K-cM)-ILx-x6ll  f  ) 

x 

and  the  residual  vectors  for  are  available  after  part  1  of  step  k  . 

The  price  paid  for  these  advantages  is  that  M  must  be  factored  al¬ 
though  it  need  not  be  invertible.  It  is  also  preferable  to  factor  K  -  OM 
if  this  is  feasible. 

II .  The  inertial  inner  product. 

Recall  that  we  want  to  use  the  operator,  or  matrix,  (K  -  OM)  1M 
because  its  eigenvalue  distribution  is  advantageous  for  both  subspace 
iteration  and  for  the  Lanczos  algorithm.  However  (K  -  OM)  -1M  is  not 
symmetric  and  it  is  often  said  that  Lanczos  requires  a  symmetric  matrix. 
That  is  not  quite  correct. 

Let  us  assume  for  the  moment  that  M  is  positive  definite.  We  will 
consider  singular  M  at  the  end  of  this  section.  The  key  observation  is 
that  (K  -  CM)”1!!  is  self-ad.1  oint  with  respect  to  the  "inertial"  inner 
product 

(u,  v)M  =  v*Mu  . 

Here  is  the  proof. 

((K  -  CM)-1Mu,  v)M  =  vV(K  -  oM)_1Mu  , 

=  v*M(K  -  CM)_1*Mu  , 

=  (u,  (X  -  OM)-1Mv)m  .  □ 


Here  is  the  algorithm. 


As  usual  lullM  =  /(u,u)M  . 


Lanczos  Algorithm  [  C  must  be  self-adjoint  with  respect  to  (*,  * )  ]  . 

Pick  r.,  ¥=  0  and  compute  g  =  Ur.,  II,.  .  Set  q_  =  0  . 

1  1  1  M  0 

For  j  =  1,  2,  ...  until  convergence  repeat 

1-  t  =  rj/ej 

2-  t 

3-  °j  -  (v  Vm 
k-  ■  uj  -  Vj 

5'  6J.l  '  ,rJ+l,M 

6.  Test  for  convergence. 

In  our  case  C  =  (K  -  crM)-1M  . 

This  formulation  is  not  widely  appreciated.  It  was  given  explicitly  in 
[van  Kats  and  van  der  Vorst,  1977].  An  alternative,  but  more  expensive 
version  of  the  algorithm  is  given  in  [Parlett,  1980,  p.  321+]. 

The  numbers  ,  8^  computed  in  the  course  of  the  algorithm  are  put 

into  a  tridiagonal  matrix 


= 


6  a 
2  2 


a 

3  3 


whose  order  j  grows  by  one  at  each  step 


The  eigenvalues  6.  ,  i  =  1,  j  of  T.  are  the  Pitz  approximations 

1  J 

to  the  eigenvalues  of  C  =  (K  -  M)  Hi  .  As  j  increases  the  outer 
values  0.  quickly  converge  to  the  outer  eigenvalues  of  (K  -  oM)  Hi 
and  these  may  be  converted,  by  (2  -  l)  ,  into  the  required  eigenvalues 
of  (K,  M)  close  to  0  . 

Let  us  consider  now  the  case  when  M  is  singular  but  K  is  positive 
definite.  It  follows  that  (K,  M)  has  one  or  more  infinite  eigenvalues 
while  C  =  (K  -  oM)  Hi  has  one  or  more  zero  eigenvalues.  Formally  the 
bilinear  expression  (•  ,  •  is  not  a  true  inner  product,  however  when 
restricted  to  the  space  corresponding  to  the  finite  eigenvalues  of  (X,  M) 
it  is  an  inner  product.  The  algorithm  given  above  automatically  keeps  the 
q^  ,  i  >  1  ,  in  this  subspace.  It  is  necessary  to  take  r  =  Mr^  to 
ensure  that  q^  is  also  in  the  subspace.  This  is  the  only  modification 
needed  when  M  is  positive  semidefinite. 

The  ghost  values  which  can  afflict  the  Ritz  approximations  for  the 
pencil  (K,  M)  cannot  affect  either  (L(K  -  cM)  1LT,  I)  or  (C,  I)  . 
Consider  again  the  example  in  Section  1.  The  Rayleigh  quotient  of  x  with 
resDect  to  both  pencils  is  p  =  H  (v~v  +  t-  1  n)  .  For  any  a  £  (x  x  ) 

2  X-a  x30-a  1  2 

we  find  P  ^  " — — )  while  the  algorithm  only  looks  at  those  Ritz 

Li-U  '•o-o 


values  outside  this  interval.  The  fact  that  Ritz  approximations  for 
inverted  operators  prevent  the  occurrence  of  ghost  vectors  was  pointed 
out  in  [Scott,  1981]. 
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h.  THE  LANTECH  ALlCP.ITnM 

Figures  1  and  2  show  a  comparison  of  the  simple  Lanczos  algorithm  and  sub¬ 
space  iteration  on  a  structural  problem.  Comparison  was  stopped  when 
subspace  iteration  exhausted  our  resources. 

To  be  fair  the  simple  Lanczos  algorithm  should  be  ccmpared  with  the  simple  ;,...er 
method  and  the  block  Lanczos  algorithm  should  be  compared  with  supspace  ir.teraticn 
using  the  same  block  size.  We  new  give  a  theoretical  comparison  to  complement  the 
the  practical  one. 

Let  A  have  eighenvalues  -V  =  i  for  i  =  0,  . ..,  n  +  1  .  Thus  A  is  (n  +  2 
by  (n  +  2)  .  Table  1  gives  the  number  steps  m  required  to  reduce  the  error  in 
the  approximation  to  zn+1  below  1"  of  its  original  value.  The  eigenvalue  error 

-  0  will  then  have  been  reduced  to  .01"  of  its  original  value.  We  assume 
that  the  initial  vector  makes  an  angle  of  1*5°  with  the  dominant  eigenvector.  The 
values  for  Lanczos  are  overestimates.  The  eigenvalue  distribution  is  a  difficult 
one  but  it  illustrates  the  power  of  the  Lanczos  algorithm. 

Table  1.  Number  of  steps  m  to  reduce  error. 


n 

Power  Method 

ljanczos  i 

2 

1CT 

163 

1 

32 

103 

16C7 

no 

iou 

16051 

380 

n 

n  £n  100 

vn  In  200 

The  superiority  of  even  the  simple  Lanczos  algorithm  over  subspace  interaticn 
can  be  explained  as  follows.  Cuppose  that  subspace  iteration  uses 
a  subspace  of  dimension  p  and  runs  for  J  steps .  It  produces 
Rayleigh  Ritz  approximation  from  a  subspace  of  dimension  p  at 
every  step.  Cn  the  other  hand,  with  the  same  computational  effort 


a 


i '  i  iiiniiitteMr'iiiM  ii  n  r  i'  •  l 


?  are  9 

and  using  the  same  inverted  operator 

Lanccos  produces,  in  exact  arithmetic,  the  Ray lei. eh  Fits  approximations 
from  a  sutspece  of  dimension  jp  .  Lar.czos  never  f  or  rets  an;/  vector  that 
has  been  computed  in  the  iteration,  although  the  vector  itself  may  have 
been  discarded  the  essential  information  is  cleverly  condensed  in  the 
tridiagonal  matrix  T  . 

There  is  no  space  here  to  give  a  detailed  account  of  the  Lanczos  algorithm 
and  the  reader  is  referred  to  [Farlett,  1980 ]  for  more  information. 

Between  1950  and  1972  the  algorithm  had  a  poor  reputation,  not  because 
it  was  poor  but  because  it  was  poorly  understood.  As  a  result  of  the 
analysis  of  Paige  the  situation  was  remedied  and  it  did  not  take  long  for 
some  good  implementations  to  appear.  One  important  point  is  that  fairly 
cheap  error  bounds  can  be  computed  during  the  iteration.  This  permits  the 
correct  information  to  be  extracted  from  the  computed  quantities  and  it 
permits  the  algorithm  to  be  stopped  as  soon  as  possible. 

Roundoff  error  has  a  significant  impact  on  the  Lancccs  algorithm.  Its 
effect  is  not  to  prevent  convergence  but  only  to  delay  it  a  little  and 
also  to  make  the  extraction  of  eigenvalues  from  the  Hits  values  somewhat  mere 
complicated. 

The  importance  of  inverting  the  proper  operator,  implicitly  of  course, 
as  discussed  in  Section  3,  is  well  illustrated  by  the  Lanczcs  process. 

There  are  many  problems  in  which  n  and  M  are  both  positive  definite,  M  is 
diagonal  or  nearly  so,  while  K  has  a  significant  bandwidth  (such  as  vn  )  . 
Most  text  books  and  references  suggest  the  following  reduction  tc  standard  form 


(K,  M)  -  (M' 


■1/2.., 


or  very 


Usually  the  eigenvalues  near  0  are  wanted  and  o  may  well  be  0 
small  relative  to  the  largest  eigenvalue.  The  eigenvalue  distribution 
typically  is 


This  is  difficult  for  Lanczos  and  disastrous  for  the  power  method  with 

the  same  operator.  The  only  attraction  of  this  approach  for  the  large  problems 

that  no  matrix  factorization  of  K  -  CM  is  needed. 

The  temptation  to  avoid  such  factoring  should  be  resisted.  Even  if 
the  cost  is  high  the  reduction  of  K  -  oM  to  LDL*  is  amply  rewarded  as  sug¬ 
gested  in  Section  3.  For  the  operator  «^2(K  -  OM)”^  M-/2  the  spectrum  is 
inverted  to: 


0  (X-cT1 


This  distribution  is  favorable  for  both  subspace  iteration  and  Lanczos  but 
Lanczos  takes  stronger  advantage  it. 

STORAGE  CONSIDERATIONS.  Lanczos  requires  in  the  fast  store,  (i)  two  or  three 
n-vectors  for  the  iteration,  (ii)  two  vectors  of  modest  length  (say  5*n  ) 

to  hold  the  tridiagonal  T  ,  (iii)  a  small  array  to  hold  Ritz  values. 

J 

It  is  also  desireable,  but  not  absolutely  necessary  to  hold  in  fast  storage 
whatever  is  necessary  to  form  Mq  and  to  solve  (K  -  oM)r  =  Mq  .  Often 
this  will  be  M  and  the  triangular  factors  L  and  D  of  K  -  cM  .  The 
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attractive  feature  of  Lanczos  is  that  each  Lanczos  vector  can  be  put  into 
secondary  storage  (a  disk  perhaps)  one  step  after  it  is  created.  These 
vectors  are  not  needed  again  to  compute  eigenvalues  but  they  are  required 
at  the  end,  after  convergence,  to  form  the  eigenvectors.  They  can  also 
be  used  from  time  to  time  to  improve  the  orthogonality  among  subsequent 
Lanczos  vectors  and  thus  to  hasten  convergence.  This  device  is  called 
selective  orthogonalizatior.  and  is  described  in  [Parlett,  19801. 

SHIFTS  OF  ORIGIN.  Sometimes  it  is  feasible  to  factor  K  -  oM  more  than 
once,  perhaps  for  five  different  values  of  a  .  This  facility  helps  both 
subspace  iteration  and  Lanczos  because  the  big  interval  [ ,  T^]  can  be 
split  into  five  smaller  ones  and  either  algorithm  can  be  used  to  find  just 
the  eigenvalues  in  the  subinterval.  The  choice  for  the  five  origin  shifts 
can  be  determined  during  the  computation.  Since  Lanczos  has  available  more 
information  on  the  spectrum  than  does  subspace  iteration  it  can  make  some¬ 
what  better  choices  for  the  o  , 

5.  SOLVING  (K  -  OM)r  =  Mu  . 

Structural  engineers  are  very  fond  of  what  they  call  (variable)  band- 
solvers  .  For  two  dimensional  problems  K  has  a  modest  bandwidth  and  L 
inherits  this  structure.  Each  partial  column  of  L  is  thought  of  as  a 
conventional  vector  from  the  top  nonzero  element  in  the  column  down  to  the 
diagonal.  In  other  words,  zero  elements  within  rhe  band  take  up  storage 
space .  The  gain  is  a  beautifully  simple  data  structure  for  L  but  for 


three  dimensional  problems  and  for  complicated  two  dimensional  structures 
the  waste  of  storage  space  is  too  expensive.  More  complicated  but  more 
effcient  data  structures  are  called  for. 

Considerable  progress  has  been  made  in  this  area  of  sparse  matrix 
technology  and  users  should  be  aware  of  the  enormous  gains  to  be  had  from 
using  an  ordering  algorithm  which  effectively  permutes  rows  and  columns  of 
K  and  M  to  nearly  minimize  the  number  of  nonzero  elements  in  L  . 
Important  algorithms  for  this  task  are  (i)  the  minimum  degree  algorithm, 
(ii)  nested  dissection,  (iii)  the  Gibbs-Poole-Stockmeyer  profile  reducer. 
Usually  these  are  all  better  than  the  standard  reverse  Cuthill-McKee 
algorithm.  It  is  probably  best  to  use  some  of  the  good  sparse  solver 
codes  now  available  rather  than  attempt  to  program  the  algorithms  one¬ 
self.  The  three  best  known  codes  are  the  Harwell  library  MAI 7  routines, 
SPARSPAK  by  Alan  George  and  J.  Liu  at  Univ.  of  Waterloo,  and  the  Yale 
Sparse  Matrix  Package.  The  most  recent  reviews  are  [Duff,  1979]  and  [Duff 
and  Stewart,  1979]. 

It  is  important  to  remember  that  the  linear  systems  to  be  solved 
should  be  indefinite  but  will  have  most  of  their  eigenvalues  positive. 

6.  DAVIDSON'S  METHOD 
/  li  . 

For  very  large  (n  >  10  )  problems  which  occur  in  atomic  and 
molecular  orbit  calculations  the  matrices  are  not  very  sparse  but  they 
are  quite  special  in  the  sense  that  the  eigenvector  matrix  is  fairly  close 
to  the  identity.  Davidson's  method  [Davidson,  1975],  which  we  have  no 
space  to  describe  here,  is  based  on  perturbation  theory  and  seems  to  be 
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very  effective.  However,  both  theoretical  and  practical  considerations 
[Kalamboukis ,  1980]  suggest  that  it  cannot  be  regarded  as  a  general  purpose 
eigenvalue  method. 
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For  all  beams  and  columns: 


Young's  Modulus 
Mass  Density 
Area 

Moment  of  Inertia 


=1.0  KN/m2 
=1.0  Kg/m3 
=  1.0m2 
=  1.0m4 


No.  of  Beam  Elements 
No.  of  Modes  =  36 
Total  Mo.  of  D.O.F. 
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Figure  lb.  Eigenvalues  of  the  above  system 
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